4.1.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- dat
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "bw_gain"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable     n  mean    sd
##   <chr>     <fct>    <dbl> <dbl> <dbl>
## 1 CTRL      bw_gain     12  17.1  2.26
## 2 PFOS      bw_gain     12  17.2  4.04
## 3 VAN       bw_gain     12  17.6  2.69
## 4 VAN+PFOS  bw_gain     12  17.5  2.65

4.1.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R01             1 no    no     310   321.  325.  339.  350    354
## 2 PFOS      R30            18 yes   no     258.  262.  265.  271.  270.   274
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.951  0.0457

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44     0.530 0.664

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers, has equal variance and is normally distributed without the outliers according to Shapiro-Wilk test. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.1.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd     F     p p<.05   ges
## 1 treatment   3  44 0.064 0.979       0.004

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term   group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
## * <chr>  <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl> <dbl> <chr>       
## 1 treat… CTRL   PFOS            0   0.111     -3.15      3.37 1     ns          
## 2 treat… CTRL   VAN             0   0.463     -2.79      3.72 0.981 ns          
## 3 treat… CTRL   VAN+P…          0   0.374     -2.88      3.63 0.99  ns          
## 4 treat… PFOS   VAN             0   0.353     -2.90      3.61 0.991 ns          
## 5 treat… PFOS   VAN+P…          0   0.264     -2.99      3.52 0.996 ns          
## 6 treat… VAN    VAN+P…          0  -0.0893    -3.35      3.17 1     ns

4.1.3.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- subset(dat, !is.na(cecum_norm))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "cecum_wt"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable     n  mean    sd
##   <chr>     <fct>    <dbl> <dbl> <dbl>
## 1 CTRL      cecum_wt    11  5.08  1.13
## 2 PFOS      cecum_wt    12  4.70  1.23
## 3 VAN       cecum_wt    12  9.41  1.27
## 4 VAN+PFOS  cecum_wt    11  9.97  1.14

4.1.3.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 6 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R01             1 no    no     310   321.  325.  339.  350    354
## 2 PFOS      R25            13 yes   no     339.  340.  353.  364.  348.   358
## 3 VAN       R15            27 no    yes    268.  277.  283.  290.  296.   300
## 4 VAN+PFOS  R43            43 yes   yes    292.  301.  300.  313.  316.   322
## 5 VAN+PFOS  R44            44 yes   yes    261.  269.  277   284.  287.   296
## 6 VAN+PFOS  R47            47 yes   yes    242.  249.  255.  263.  267.   271
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains six not critical outliers.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.954  0.0655

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    42    0.0879 0.966

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that normalised cecum weight data has six non-critical outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.1.3.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(cecum_wt ~ pfos*van)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##     Effect DFn DFd       F        p p<.05   ges
## 1     pfos   1  42   0.066 7.99e-01       0.002
## 2      van   1  42 185.208 5.41e-17     * 0.815
## 3 pfos:van   1  42   1.730 1.96e-01       0.040

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term      group1 group2   null.value estimate conf.low conf.high    p.adj
## * <chr>     <chr>  <chr>         <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
## 1 treatment CTRL   PFOS              0   -0.374   -1.71      0.961 8.77e- 1
## 2 treatment CTRL   VAN               0    4.34     3.00      5.67  3.69e-10
## 3 treatment CTRL   VAN+PFOS          0    4.89     3.53      6.26  2.36e-11
## 4 treatment PFOS   VAN               0    4.71     3.41      6.02  2   e-11
## 5 treatment PFOS   VAN+PFOS          0    5.27     3.93      6.60  2.39e-12
## 6 treatment VAN    VAN+PFOS          0    0.554   -0.780     1.89  6.85e- 1
## # ℹ 1 more variable: p.adj.signif <chr>

4.2.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- subset(dat, !is.na(cecum_norm))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "cecum_norm"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable       n  mean    sd
##   <chr>     <fct>      <dbl> <dbl> <dbl>
## 1 CTRL      cecum_norm    11 1     0.148
## 2 PFOS      cecum_norm    12 0.944 0.164
## 3 VAN       cecum_norm    12 1.88  0.24 
## 4 VAN+PFOS  cecum_norm    11 2.07  0.201

4.2.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 VAN       R15            27 no    yes    268.  277.  283.  290.  296.   300
## 2 VAN       R24            36 no    yes    281.  286.  294.  305.  309.   312
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.984   0.753

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    42     0.416 0.742

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that normalised cecum weight data has two outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.2.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd       F        p p<.05   ges
## 1 treatment   3  42 106.226 1.21e-19     * 0.884

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term      group1 group2   null.value estimate conf.low conf.high    p.adj
## * <chr>     <chr>  <chr>         <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
## 1 treatment CTRL   PFOS              0  -0.0562  -0.271      0.158 8.96e- 1
## 2 treatment CTRL   VAN               0   0.884    0.669      1.10  1.42e-12
## 3 treatment CTRL   VAN+PFOS          0   1.07     0.850      1.29  1.06e-12
## 4 treatment PFOS   VAN               0   0.940    0.730      1.15  1.09e-12
## 5 treatment PFOS   VAN+PFOS          0   1.13     0.911      1.34  1.06e-12
## 6 treatment VAN    VAN+PFOS          0   0.186   -0.0287     0.400 1.1 e- 1
## # ℹ 1 more variable: p.adj.signif <chr>

4.3.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove NA in the data column
dat.clean <- subset(dat, !is.na(liver_norm))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "liver_wt"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable     n  mean    sd
##   <chr>     <fct>    <dbl> <dbl> <dbl>
## 1 CTRL      liver_wt    12  11.2 1.12 
## 2 PFOS      liver_wt    12  12.0 1.51 
## 3 VAN       liver_wt    12  10.5 0.973
## 4 VAN+PFOS  liver_wt    12  11.3 1.47

4.3.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 1 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R01             1 no    no      310  321.  325.  339.   350   354
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains one non-critical outlier.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.989   0.922

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44      1.50 0.227

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that normalised liver weight data has one non-critical outlier, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.3.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(liver_wt ~ pfos*van) #FORMULA
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##     Effect DFn DFd        F     p p<.05      ges
## 1     pfos   1  44 4.977000 0.031     * 1.02e-01
## 2      van   1  44 3.492000 0.068       7.40e-02
## 3 pfos:van   1  44 0.000905 0.976       2.06e-05

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA) #FORMULA
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term  group1 group2 null.value estimate conf.low conf.high  p.adj p.adj.signif
## * <chr> <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl>  <dbl> <chr>       
## 1 trea… CTRL   PFOS            0    0.820   -0.587     2.23  0.414  ns          
## 2 trea… CTRL   VAN             0   -0.708   -2.11      0.699 0.541  ns          
## 3 trea… CTRL   VAN+P…          0    0.135   -1.27      1.54  0.994  ns          
## 4 trea… PFOS   VAN             0   -1.53    -2.93     -0.121 0.0287 *           
## 5 trea… PFOS   VAN+P…          0   -0.685   -2.09      0.722 0.568  ns          
## 6 trea… VAN    VAN+P…          0    0.843   -0.565     2.25  0.39   ns

4.3.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "µg/g") + #,limits = c(8,17),breaks = seq(8,17,2)
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE) #, y.position = c(1.35,1.4,1.45,1.5)
p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Liver weights at day 8

4.4.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove NA in the data column
dat.clean <- subset(dat, !is.na(liver_norm))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "liver_norm"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable       n  mean    sd
##   <chr>     <fct>      <dbl> <dbl> <dbl>
## 1 CTRL      liver_norm    12 1     0.055
## 2 PFOS      liver_norm    12 1.08  0.055
## 3 VAN       liver_norm    12 0.934 0.044
## 4 VAN+PFOS  liver_norm    12 1.03  0.065

4.4.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Data contains zero outliers.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.985   0.778

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44     0.430 0.733

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that normalised liver weight data has two outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.4.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd      F        p p<.05  ges
## 1 treatment   3  44 15.909 3.79e-07     * 0.52

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term      group1 group2   null.value estimate conf.low conf.high       p.adj
## * <chr>     <chr>  <chr>         <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
## 1 treatment CTRL   PFOS              0   0.0846   0.0246   0.145   0.0027     
## 2 treatment CTRL   VAN               0  -0.0663  -0.126   -0.00623 0.0254     
## 3 treatment CTRL   VAN+PFOS          0   0.0352  -0.0249   0.0953  0.409      
## 4 treatment PFOS   VAN               0  -0.151   -0.211   -0.0909  0.000000182
## 5 treatment PFOS   VAN+PFOS          0  -0.0494  -0.110    0.0107  0.14       
## 6 treatment VAN    VAN+PFOS          0   0.102    0.0415   0.162   0.000269   
## # ℹ 1 more variable: p.adj.signif <chr>

4.4.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% difference",limits = c(0.75,1.5),breaks = seq(0.75,1.5,0.25), labels = function(x) paste0(x*100, "%")) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE, y.position = c(1.35,1.4,1.45,1.5))
p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Normalized liver weights at day 8

5.1.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove rows with NA
dat.clean <- subset(dat, !is.na(pfos_serum4_ugml))
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "pfos_serum4_ugml"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable             n   mean    sd
##   <chr>     <fct>            <dbl>  <dbl> <dbl>
## 1 CTRL      pfos_serum4_ugml    11  0     0.001
## 2 PFOS      pfos_serum4_ugml    12  9.17  2.01 
## 3 VAN       pfos_serum4_ugml    11  0.001 0.001
## 4 VAN+PFOS  pfos_serum4_ugml    10 10.0   1.18

5.1.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 5 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R06             6 no    no     274.  276.  285.  296.  297.   302
## 2 CTRL      R11            11 no    no     271.  278.  287.  294.  293.   299
## 3 PFOS      R29            17 yes   no     239.  243.  248.  255.  259.   267
## 4 VAN       R13            25 no    yes    218.  222.  228.  234   231.   241
## 5 VAN       R21            33 no    yes    262.  268.  274.  285.  281.   291
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic    p.value
##   <chr>                <dbl>      <dbl>
## 1 residuals(model)     0.807 0.00000412

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1     3    40      7.48 0.000436

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers and has equal variance, however falls short on the Shapiro-Wilk test of normality and is therefore not normally distributed. Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.1.1.3 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.                  n statistic    df           p method        
## * <chr>            <int>     <dbl> <int>       <dbl> <chr>         
## 1 pfos_serum4_ugml    44      35.4     3 0.000000101 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.                  n effsize method  magnitude
## * <chr>            <int>   <dbl> <chr>   <ord>    
## 1 pfos_serum4_ugml    44   0.809 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc

## # A tibble: 6 × 9
##   .y.           group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
## * <chr>         <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 pfos_serum4_… CTRL   PFOS      11    12    3.86   1.15e-4 1.85e-4 ***         
## 2 pfos_serum4_… CTRL   VAN       11    11    0.0172 9.86e-1 9.86e-1 ns          
## 3 pfos_serum4_… CTRL   VAN+P…    11    10    4.53   5.87e-6 1.91e-5 ****        
## 4 pfos_serum4_… PFOS   VAN       12    11   -3.84   1.23e-4 1.85e-4 ***         
## 5 pfos_serum4_… PFOS   VAN+P…    12    10    0.863  3.88e-1 4.66e-1 ns          
## 6 pfos_serum4_… VAN    VAN+P…    11    10    4.51   6.36e-6 1.91e-5 ****

5.1.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "ug/mL",limits = c(0,20),breaks = seq(0,20,5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE, y.position = c(14,17,15,14))
p

## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 16 rows containing missing values (`geom_point()`).

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
## Removed 16 rows containing missing values (`geom_point()`).

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
## Removed 16 rows containing missing values (`geom_point()`).

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS conc. in ug/mL at day 4

5.1.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_serum4_mg"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes")

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_serum4_mg))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.1.2.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable       statistic     p
##   <chr>     <chr>              <dbl> <dbl>
## 1 PFOS      pfos_serum4_mg     0.955 0.717
## 2 VAN+PFOS  pfos_serum4_mg     0.930 0.446

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    20    0.0274 0.870

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

No outliers were identified. Data is normally distributed and has equal variance. Hence we use t-test.

5.1.2.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1  -0.0171     0.165     0.183 pfos_s… PFOS   VAN+P…    12    10     -1.23 0.232
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

The output provides:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• statistic: Test statistic used to compute the p-value.

• df: degrees of freedom.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size. It corresponds to the estimated mean or difference in means depending on whether it was a one-sample test or a two-sample test.

• estimate1, estimate2: show the mean values of the two groups, respectively, for independent samples t-tests.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval.

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.            group1 group2   effsize    n1    n2 magnitude
## * <chr>          <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum4_mg PFOS   VAN+PFOS  -0.528    12    10 moderate

5.1.2.4 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "mg PFOS",limits = c(0,0.30),breaks = seq(0,0.30,0.1)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE, y.position = c(0.28))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Total mg PFOS in blood volume

5.1.3.1 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 1 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 PFOS      R29            17 yes   no     239.  243.  248.  255.  259.   267
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable        statistic     p
##   <chr>     <chr>               <dbl> <dbl>
## 1 PFOS      pfos_serum4_pct     0.937 0.456
## 2 VAN+PFOS  pfos_serum4_pct     0.939 0.547

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    20      1.43 0.246

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

5.1.3.2 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1   -0.626      6.75      7.37 pfos_s… PFOS   VAN+P…    12    10     -1.17 0.254
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

The output provides:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• statistic: Test statistic used to compute the p-value.

• df: degrees of freedom.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size. It corresponds to the estimated mean or difference in means depending on whether it was a one-sample test or a two-sample test.

• estimate1, estimate2: show the mean values of the two groups, respectively, for independent samples t-tests.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval.

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.             group1 group2   effsize    n1    n2 magnitude
## * <chr>           <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum4_pct PFOS   VAN+PFOS  -0.503    12    10 moderate

5.1.3.3 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% of total dosed PFOS", limits = c(3,10),breaks = seq(3,10,1)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE) #, y.position = c(1.35,1.4,1.45,1.5))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS serum day 4 in pct. of total

5.2.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove rows with NA
dat.clean <- subset(dat, !is.na(pfos_serum8_ugml))
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "pfos_serum8_ugml"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable             n   mean     sd
##   <chr>     <fct>            <dbl>  <dbl>  <dbl>
## 1 CTRL      pfos_serum8_ugml    12  0.016  0.03 
## 2 PFOS      pfos_serum8_ugml    12 36.3   15.5  
## 3 VAN       pfos_serum8_ugml    12  0.011  0.021
## 4 VAN+PFOS  pfos_serum8_ugml    12 32.2   10.7

5.2.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 7 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R05             5 no    no     214   219.  222.  229.  231.   237
## 2 CTRL      R09             9 no    no     273.  278.  290.  290.  296.   302
## 3 CTRL      R12            12 no    no     297.  304.  316.  322.  324.   334
## 4 VAN       R14            26 no    yes    246.  256.  260.  267.  270.   274
## 5 VAN       R16            28 no    yes    256.  260   268.  275.  273    279
## 6 VAN       R24            36 no    yes    281.  286.  294.  305.  309.   312
## 7 VAN+PFOS  R47            47 yes   yes    242.  249.  255.  263.  267.   271
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic    p.value
##   <chr>                <dbl>      <dbl>
## 1 residuals(model)     0.823 0.00000461

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic         p
##   <int> <int>     <dbl>     <dbl>
## 1     3    44      8.92 0.0000993

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers and has equal variance, however falls short on the Shapiro-Wilk test of normality and is therefore not normally distributed. Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.2.1.3 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.                  n statistic    df            p method        
## * <chr>            <int>     <dbl> <int>        <dbl> <chr>         
## 1 pfos_serum8_ugml    48      37.3     3 0.0000000402 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.                  n effsize method  magnitude
## * <chr>            <int>   <dbl> <chr>   <ord>    
## 1 pfos_serum8_ugml    48   0.779 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc

## # A tibble: 6 × 9
##   .y.           group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
## * <chr>         <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 pfos_serum8_… CTRL   PFOS      12    12    4.37   1.26e-5 3.78e-5 ****        
## 2 pfos_serum8_… CTRL   VAN       12    12   -0.0899 9.28e-1 9.28e-1 ns          
## 3 pfos_serum8_… CTRL   VAN+P…    12    12    4.17   3.02e-5 4.53e-5 ****        
## 4 pfos_serum8_… PFOS   VAN       12    12   -4.46   8.32e-6 3.78e-5 ****        
## 5 pfos_serum8_… PFOS   VAN+P…    12    12   -0.195  8.46e-1 9.28e-1 ns          
## 6 pfos_serum8_… VAN    VAN+P…    12    12    4.26   2.03e-5 4.05e-5 ****

5.2.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "ug/mL",limits = c(0,80),breaks = seq(0,80,10)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE, y.position = c(72,80,75,72))
p

## Warning: Removed 11 rows containing missing values (`geom_point()`).

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 11 rows containing missing values (`geom_point()`).

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 11 rows containing missing values (`geom_point()`).

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS conc. in ug/mL at day 8

5.2.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_serum8_mg"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes")

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_serum8_mg))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.2.2.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable       statistic     p
##   <chr>     <chr>              <dbl> <dbl>
## 1 PFOS      pfos_serum8_mg     0.890 0.119
## 2 VAN+PFOS  pfos_serum8_mg     0.902 0.168

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    22      1.98 0.173

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

No outliers were identified. Data is normally distributed and has equal variance. Hence we use t-test.

5.2.2.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1   0.0931     0.717     0.624 pfos_s… PFOS   VAN+P…    12    12     0.853 0.403
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.            group1 group2   effsize    n1    n2 magnitude
## * <chr>          <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum8_mg PFOS   VAN+PFOS   0.348    12    12 small

5.2.2.4 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "mg PFOS",limits = c(0,2),breaks = seq(0,2,0.5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE, y.position = c(1.75))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Total mg PFOS in blood volume

5.2.3.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_serum8_pct"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes" & !rat_name == "R47")

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_serum8_pct))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.2.3.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable        statistic      p
##   <chr>     <chr>               <dbl>  <dbl>
## 1 PFOS      pfos_serum8_pct     0.867 0.0594
## 2 VAN+PFOS  pfos_serum8_pct     0.899 0.181

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    21      2.86 0.106

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

No outliers were identified. Data is normally distributed and has equal variance. Hence we use t-test.

5.2.3.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1     2.05      11.8      9.74 pfos_s… PFOS   VAN+P…    12    11      1.29 0.212
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.             group1 group2   effsize    n1    n2 magnitude
## * <chr>           <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum8_pct PFOS   VAN+PFOS   0.537    12    11 moderate

5.2.3.4 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% of total dosed PFOS", limits = c(5,25),breaks = seq(5,25,5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE, y.position = c(24))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS serum day 8 in pct. of total

5.3.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_change48_pct"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes") # add following to subset() to remove the outliers: & !rat_name %in% c("R47","R27"))

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_change48_pct))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.3.1.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 PFOS      R27            15 yes   no     270.  280.  284.  291.  290.   296
## 2 VAN+PFOS  R47            47 yes   yes    242.  249.  255.  263.  267.   271
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Any extreme outliers can be bad samples or errors in data entry. If outliers, compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable          statistic      p
##   <chr>     <chr>                 <dbl>  <dbl>
## 1 PFOS      pfos_change48_pct     0.894 0.132 
## 2 VAN+PFOS  pfos_change48_pct     0.805 0.0167

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    20     0.639 0.434

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

Two outliers were identified (sample for R27 and R47). Analysis result and test method is similar with and without outliers. Data is normally distributed and has equal variance. Hence we use t-test.

5.3.1.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1     85.1      341.      256. pfos_c… PFOS   VAN+P…    12    10      1.11 0.281
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.               group1 group2   effsize    n1    n2 magnitude
## * <chr>             <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_change48_pct PFOS   VAN+PFOS   0.475    12    10 small

5.3.1.4 Conclusion

5.3.1.5 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create point plot with mean and SD
data_summary <- function(x) {
  m <- mean(x)
  ymin <- m-sd(x)
  ymax <- m+sd(x)
  return(c(y=m,ymin=ymin,ymax=ymax))
}
data_summary_collapsed <- function(x) {
  m <- mean(x)
  ymin <- m
  ymax <- m
  return(c(y=m,ymin=ymin,ymax=ymax))
}

p <- ggplot(dat.clean, aes(x = .data[[PREDICTOR]], y = .data[[OUTCOME]], color = .data[[PREDICTOR]])) +
  stat_summary(fun.data = data_summary_collapsed, geom = "crossbar", color = "black", width = 0.5, linewidth = 0.3) +
  stat_summary(fun.data = data_summary, geom = "errorbar", color = "black", width = 0.15, linewidth = 0.5) +
  geom_point(position = position_jitterdodge(dodge.width = 0.6, jitter.width = 0.4), size = 2, colour = "black", shape = 21, stroke = 0.5, aes(fill = treatment)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% change", limits = c(100,900),breaks = seq(100,900,100), labels = function(x) paste0(x, "%")) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment") +
  theme_pubr()
p

# Alternative: Create boxplot
# p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
#           fill = PREDICTOR,
#           add =  "jitter",
#           add.params = list(size = 1)) +
#   scale_fill_manual(values = params$COL) +
#   scale_y_continuous(name = "% change", limits = c(100,900),breaks = seq(100,900,100)) +
#   labs(fill = "Treatment") +
#   scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE) #, y.position = c(1.35,1.4,1.45,1.5))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

p

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 70, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS serum level change from day 4 to 8 in pct.

5.3.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Color scheme
COL <- c("#61d46b","#ffe900","#31b44b","#efc000")

# Subset data
dat.sub <- subset(dat, pfos == "yes")

# Create data frame for data representation
dat.clean <- dat.sub %>% select(rat_name, treatment, pfos_serum4_ugml, pfos_serum8_ugml) %>%
  pivot_longer(., cols = c(pfos_serum4_ugml, pfos_serum8_ugml), names_to = "data_group", values_to = "conc")

# Create column for day of sampling
dat.clean <- transform(dat.clean, "day" = ifelse(dat.clean$data_group == "pfos_serum8_ugml","d8","d4"))

# Create ID column for easier handling
for (i in dat.sub$rat_name) {
  dat.clean$ID <- paste(dat.clean$day,"_",dat.clean$treatment)
}

# Order dataframe for analysis
dat.clean <- dat.clean[order(dat.clean$day),]

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(conc))
dat.clean

##    rat_name treatment       data_group  conc day            ID
## 1       R25      PFOS pfos_serum4_ugml  8.56  d4     d4 _ PFOS
## 3       R26      PFOS pfos_serum4_ugml  9.60  d4     d4 _ PFOS
## 5       R27      PFOS pfos_serum4_ugml  7.92  d4     d4 _ PFOS
## 7       R28      PFOS pfos_serum4_ugml  8.64  d4     d4 _ PFOS
## 9       R29      PFOS pfos_serum4_ugml 12.96  d4     d4 _ PFOS
## 11      R30      PFOS pfos_serum4_ugml  8.68  d4     d4 _ PFOS
## 13      R31      PFOS pfos_serum4_ugml  5.72  d4     d4 _ PFOS
## 15      R32      PFOS pfos_serum4_ugml  7.56  d4     d4 _ PFOS
## 17      R33      PFOS pfos_serum4_ugml  8.24  d4     d4 _ PFOS
## 19      R34      PFOS pfos_serum4_ugml 11.36  d4     d4 _ PFOS
## 21      R35      PFOS pfos_serum4_ugml  9.00  d4     d4 _ PFOS
## 23      R36      PFOS pfos_serum4_ugml 11.84  d4     d4 _ PFOS
## 25      R37  VAN+PFOS pfos_serum4_ugml  9.36  d4 d4 _ VAN+PFOS
## 27      R38  VAN+PFOS pfos_serum4_ugml 11.76  d4 d4 _ VAN+PFOS
## 29      R39  VAN+PFOS pfos_serum4_ugml 10.64  d4 d4 _ VAN+PFOS
## 31      R40  VAN+PFOS pfos_serum4_ugml 12.12  d4 d4 _ VAN+PFOS
## 33      R41  VAN+PFOS pfos_serum4_ugml  9.12  d4 d4 _ VAN+PFOS
## 35      R42  VAN+PFOS pfos_serum4_ugml  9.80  d4 d4 _ VAN+PFOS
## 37      R43  VAN+PFOS pfos_serum4_ugml 10.08  d4 d4 _ VAN+PFOS
## 43      R46  VAN+PFOS pfos_serum4_ugml  9.88  d4 d4 _ VAN+PFOS
## 45      R47  VAN+PFOS pfos_serum4_ugml  8.80  d4 d4 _ VAN+PFOS
## 47      R48  VAN+PFOS pfos_serum4_ugml  8.60  d4 d4 _ VAN+PFOS
## 2       R25      PFOS pfos_serum8_ugml 43.92  d8     d8 _ PFOS
## 4       R26      PFOS pfos_serum8_ugml 45.60  d8     d8 _ PFOS
## 6       R27      PFOS pfos_serum8_ugml 69.92  d8     d8 _ PFOS
## 8       R28      PFOS pfos_serum8_ugml 21.92  d8     d8 _ PFOS
## 10      R29      PFOS pfos_serum8_ugml 26.56  d8     d8 _ PFOS
## 12      R30      PFOS pfos_serum8_ugml 47.84  d8     d8 _ PFOS
## 14      R31      PFOS pfos_serum8_ugml 21.84  d8     d8 _ PFOS
## 16      R32      PFOS pfos_serum8_ugml 29.36  d8     d8 _ PFOS
## 18      R33      PFOS pfos_serum8_ugml 22.08  d8     d8 _ PFOS
## 20      R34      PFOS pfos_serum8_ugml 52.48  d8     d8 _ PFOS
## 22      R35      PFOS pfos_serum8_ugml 29.84  d8     d8 _ PFOS
## 24      R36      PFOS pfos_serum8_ugml 23.76  d8     d8 _ PFOS
## 26      R37  VAN+PFOS pfos_serum8_ugml 36.08  d8 d8 _ VAN+PFOS
## 28      R38  VAN+PFOS pfos_serum8_ugml 31.28  d8 d8 _ VAN+PFOS
## 30      R39  VAN+PFOS pfos_serum8_ugml 25.92  d8 d8 _ VAN+PFOS
## 32      R40  VAN+PFOS pfos_serum8_ugml 23.92  d8 d8 _ VAN+PFOS
## 34      R41  VAN+PFOS pfos_serum8_ugml 21.68  d8 d8 _ VAN+PFOS
## 36      R42  VAN+PFOS pfos_serum8_ugml 40.96  d8 d8 _ VAN+PFOS
## 38      R43  VAN+PFOS pfos_serum8_ugml 25.60  d8 d8 _ VAN+PFOS
## 40      R44  VAN+PFOS pfos_serum8_ugml 37.44  d8 d8 _ VAN+PFOS
## 42      R45  VAN+PFOS pfos_serum8_ugml 25.36  d8 d8 _ VAN+PFOS
## 44      R46  VAN+PFOS pfos_serum8_ugml 34.72  d8 d8 _ VAN+PFOS
## 46      R47  VAN+PFOS pfos_serum8_ugml 59.52  d8 d8 _ VAN+PFOS
## 48      R48  VAN+PFOS pfos_serum8_ugml 23.76  d8 d8 _ VAN+PFOS

# Set names of variables
PREDICTOR <- "ID"
OUTCOME <- "conc"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   ID            variable     n  mean    sd
##   <chr>         <fct>    <dbl> <dbl> <dbl>
## 1 d4 _ PFOS     conc        12  9.17  2.01
## 2 d4 _ VAN+PFOS conc        10 10.0   1.18
## 3 d8 _ PFOS     conc        12 36.3  15.5 
## 4 d8 _ VAN+PFOS conc        12 32.2  10.7

5.3.2.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = COL)
bxp

5.3.2.3 Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 8
##   ID            rat_name treatment data_group   conc day   is.outlier is.extreme
##   <chr>         <chr>    <chr>     <chr>       <dbl> <chr> <lgl>      <lgl>     
## 1 d4 _ PFOS     R29      PFOS      pfos_serum…  13.0 d4    TRUE       FALSE     
## 2 d8 _ VAN+PFOS R47      VAN+PFOS  pfos_serum…  59.5 d8    TRUE       FALSE

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.875 0.000151

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic       p
##   <int> <int>     <dbl>   <dbl>
## 1     3    42      6.46 0.00108

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that PFOS concentrations from day 4 and 8 has two outliers, has unequal variance, and falls short on the Shapiro-Wilk test of normality (not normally distributed). Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.3.2.4 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.       n statistic    df           p method        
## * <chr> <int>     <dbl> <int>       <dbl> <chr>         
## 1 conc     46      34.4     3 0.000000165 Kruskal-Wallis